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ABSTRACT 

The ONC appears to be unusual on two grounds: The observed constellation of the 
OB-stars of the entire Orion Nebula cluster and its Trapezium at its centre implies 
a time-scale problem given the age of the Trapezium, and an IMF problem for the 
whole OB-star population in the ONC. Given the estimated crossing time of the 
Trapezium, it ought to have totally dynamically decayed by now. Furthermore, by 
combining the lower limit of the ONC mass with a standard IMF it emerges that 
the ONC should have formed at least about 40 stars heavier than 5 M Q while only 
ten are observed. Using iV-body experiments we (i) confirm the expected instability 
of the trapezium and (ii) show that beginning with a compact OB-star configuration 
of about 40 stars the number of observed OB stars after 1 Myr within 1 pc radius 
and a compact trapezium configuration can both be reproduced. These two empirical 
constraints thus support our estimate of 40 initial OB stars in the cluster. Interestingly 
a more-evolved version of the ONC resembles the Upper Scorpius OB association. The 
JV-body experiments are performed with the new C-code CATENA by integrating the 
equations of motion using the chain-multiple-regularization method. In addition we 
present a new numerical formulation of the initial mass function. 

Key words: methods: TV-body simulations - open clusters and association: individual: 
ONC - stars: kinematics - stars: mass function 



1 INTRODUCTION 

Of all O-stars 46 per cent and of all B-stars 4 per 
cent are runaways exceeding 30 km/s iStondli99f[) . Fur- 
thermore the binary fraction amo ng runaway O-stars is 
around 10 % iGies fc Boltonl Il986ft while it is more than 
50 % in young star clusters jGoodwin et al.ll200rl . This 
suggests that binaries are involved in close dynamical 
encounters leading to stellar ejections while the binary 
fraction among t he ej ected stars is decreased. Indeed, 
IClarke fc Pringiel lll992t) deduced using an analytical ap- 
proach that massive stars must form in compact small- 
N groups. The decay of non-hierarchical 3,4,5-body sys- 
tems with equal mas ses as well as a mass sp ectrum has 
been investigated by ISterzik fc Durisenl £l998). They de- 
termined the spectrum of the remnant decay products but 
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not the phase-s pace behaviour of compact few-body sys - 
tems with time. [Hoogerwer^ de Bruiine fc de Zeeuwl (|200(J) 
and iHoogerweTfTtoall (I2OO1I) were able to trace back 
the trajectories of some runaways to nearby associations. 
iRamspeck. Heber fc Moehlerl teOOll) determined the age and 
the calculated time-of-flight of early-type stars at high galac- 
tic latitudes and concluded that they can have their origin 
in the galactic disk. In the case of the runaways AE Auri- 
gae and fx Columbae, which have spatial velocities greater 
than 100 km/s, in combination with the bina ry 1 Orionis, 
iGualandris. Portegies Zwart fc Eggletonl <l2004l) have shown 
that the encounter of two binaries with high eccentricities 
2.5 Myr ago and in the co-moving vicinity of the current 
Orion Nebula Cluster (ONC) can reproduce the spatial con- 
figuration observed today. The spatial distribution of field 
OB-stars can thus be understood qualitatively using the- 
oretical stellar-dynamical methods. But to obtain a more 
complete picture we need to study the details and frequency 
of occurrence of energetic ejections from the acceleration 
centres, namely the inner regions of young star clusters. 
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Table 1. Setup data for all three models 



model 




stars 


-Mtot.OB 
M 


_£3D 

km/s 


Kyr 


runs 


4-body 


e 1 


in Tab.0 


88.4 


3.9 


12.6 


1000 


10-body 


all 


in Tab. 


167.2 


5.4 


9.2 


1000 


40-body 


4 


X Tab.Hl 


668.8 


10.7 


4.7 


1000 



Specification of the N-body systems. Mtot is the total mass in 
Mq. The velocity dispersion cr is calculated from the total mass 
placed within the initial radius of 0.025 pc corresponding to the 
ONC-TS size if virial equilibrium is assumed. t cr is the crossing 
time, and runs are the total number of integrated configurations. 



2 MOTIVATING PROBLEMS 

In the case of the Orion Nebula cluster two main discrepan- 
cies concerning the properties of its OB-stars are found: 

2.1 Existence of the Trapezium system 

Given the total mass of all four Trapezium stars 
jHillenbrandlll997l) of 88.4 Mq and their occupying space of 
nearly 0.05 pc in diameter the corresponding crossing time 
can be estimated, 

*- = Vef • (1) 

to be about 13 Kyr (Tab. 0, if the Trapezium is as- 
sumed to be a compact virialized subsystem of the ONC. 
ISterzik fc Durisenl {l998) noted that most systems in their 
decay analysis decay within dozens of crossing times. So the 
ONC-TS is expecte d to have tota lly decayed by now, if its 
age is about 1 Myr jKroupalbOol) . 

2.2 The number of OB stars 

The virial mass of the ONC is measured to be nearly 
4500 M but only about 1800 M is visible in stel- 
lar mass (Hil lenbrand fc Hartmannl [598) while cluster- 
formation models that match the ONC suggest that it may 
have formed with 10 4 stars plus brown dwarfs and that it 
is expanding now resembling a Pleiades type cluster em- 
bedded in an e xpanding association to a remarkable degree 
after 100 Myr (iKroupa et alJl2001l) . The observed mass of 
all stars heavier than 5 Mq is 167 Mq. If the canonical IMF 
(App. [bJ is normalized such that 1633 Mq are contained 
in the mass interval ranging from 0.01 up to 5 Mq and for 
three different physically possibl e upper stellar mass lim- 
its, m max », of 80, 150 Mq or +0 llWeidner fc; Kroupall2004t 
lOev fc Clarke! 120051: lFigerll2005l : lKoenll200rj) . and different 
IMF-slopes above 1 Mq the corresponding maximum stellar 
mass m max and the expected number of OB-stars formed in 
the ONC can be cal culated (Tab. El - The calculations are 
based on the IMF bv lKroupal (l200lT) but with a numerically 
more convenient description 1 fApp. lAT)l . 

Given the values in Tab.|5|the ONC should have formed 
about 38 OB-stars assuming the IMF to be canonical (03 = 

1 A utility-IMF package and CATENA including a full doc- 
umentation can be downloaded from the AlfA-webpage: 
7//www.astro.uni-bonn.de/| 



2.35). But in a thorough survey of the ONC iHillenbrandl 
(1997) lists only 10 stars weighting more than 5 Mq (Tab. 
13 . 9 of these 10 OB stars are located within a projected 
sphere of 1 pc around the Trapezium system. The remaining 
B star (5.7 Mq) is placed approximately 2.3 pc away from 
the Trapezium. 7 OB stars, including the three most massive 
stars, are located within 0.5 pc around the Trapezium in 
projection. If the IMF is steepened above 1 Mq to Q3 = 2.7 
the number of expected OB-stars decreases down to 18. But 
the expected maximum stellar mass also decreases down to 
m max = 28 Mq , whereas two observed stars are heavier. The 
existence of these stars suggests that the IMF was indeed 
normal. Note that the time-scale problem would persist even 
if we allow Q3 = 2.7. Because the IMF seems to be universal 
l|Kroupal f20021 a significant deviation from the calculated 
number of 38 stars using the canonical IMF should not be 
expected. 

As a check the number of stars heavier than 5 Mq 
can also be estimated by normalising the canonical IMF 
to the number of stars in the ma ss range 1-2 Mq in the 
cluster. Using the stellar sample of IHillenbrandl (|l997T ). the 
number of stars heavier than 5 Mq can be derived from 
the number of stars between 1 and 2 Mq (70) and noting 
that in this mass regime only the non-embedded sources 
are listed. These a mount to approximately half of all stars 
dHillenbrandlll997l) . Thus, 26 OB-stars are expected to have 
formed in the ONC. The total mass derived from this mass 
regime is 1404 M , 22 % less than the total estimated mass 
used above. Given this uncertainty (13-38 stars heavier than 
5 M ), we perform computations with 10 and 40 stars. As 
will become apparent below, 40 OB stars are our preferred 
value. 

Furthermore, if stars are drawn randomly from a uni- 
versal IMF, the number of stars heavier than 5 Mq may not 
be the expectation value of 38. The number can be smaller. 
To estimate the probability that less than k of n stars have 
masses less than 5 Mq, drawing stars from an IMF has to 
be interpreted as a Bernoulli experiment: For the mass of 
the ONC, A/onc, the total number of stars, n to t, and the 
number of stars, n>5, heavier than 5 Mq can be calculated. 
If one star is drawn from the IMF, the probability to get a 
star heavier than 5 Mq is 

P = n >5 /n to t ■ (2) 

This experiment is repeated ntot times. So the probability 
to have a star heavier than 5 Mq k times is given by the 
Bernoulli-distribution, 

P{k) = (^P'il-Pr ^ ■ (3) 

Because the event probability is small and the number of 
experiments large the Poissonian limit can be applied. The 
probability is approximately 

p(fc) = f^e-" , (4) 
fi. 

where u = p n to t = n> 5 . So the total probability to get k or 
fewer OB-stars is 

P(<fc)=$>(0 ■ (5) 

1 = 

The probability to get 20, 10 or fewer OB-stars for two differ- 
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Table 2. The number of expected OB stars and maximum stellar 
mass in the Orion Nebula cluster. 



03 


m ma »/Mg : 


+0O 


150 


80 


obs. 


2.35 


m max /M0 


76.4 


59.5 


46.8 


45.7 


2.35 


N >5 


38.7 


38.3 


37.8 


10 


2.35 


M tot /M Q 


2103 


2076 


2048 


1800 


2.7 


m max /M 


28.6 


27.7 


26.1 


45.7 


2.7 


N >5 


18.4 


18.4 


18.3 


10 


2.7 


A/tot /M Q 


1799 


1797 


1794 


1800 



Observed iHillenbrandll997t> and expected maximum stellar mass 
("imax), number of stars more massive than 5 Mq (Af>s), and 
total initial mass (Mtot) for the Orion Nebula Cluster in depen- 
dence of three different physically possible upper stellar mass lim- 
its, m max *, and two different IMF-slopes, 03, for the mass range 
from 1 Mq up to m ma x- The mass rang e less than 1 Mq is de- 
scribed by a Kroupa-IMF iKrouDall200ll) . 



Table 3. Probability of a deviation from a canonical IMF 



M O nc/M 


1800 


1800 


2200 


2200 


«lmax*/M0 


80 


150 


80 


150 


"tot 


5209 


5144 


6337 


6251 




33 


33 


41 


41 


P(k < 10) 


2.8 ■ 10~ 6 


2.8 • 10~ 6 


7.6 ■ 10~ 9 


7.6 • 10" : 


P(k < 20) 


1.0 ■ 10~ 2 


1.0 • 10~ 2 


2.2 ■ 10~ 4 


2.2 • 10~ 



3 INTEGRATOR 

To investigate the dynamics of the OB-stars in the ONC 
we perform direct N-body integrations. Because close en- 
counters with high eccentricity are very frequent in com- 
pact few-body systems due to the grainy potential, a multi- 
ple regularization technique is required to reduce energy er- 
rors and speed-up the calculations. We combined in our own 
code (catena ) the very efficient CHAlN- r egula r izatio n for- 
malism developed bv lMikkola fc Aarsethl (Il99fi ll99oT) with 
an embedded Runge-Kutta method of 8(9)-th order u sing 
a coefficient-set published bv lPrince fc Dormandl dl98ll) . in- 
stead of the Aarseth-CHAIN-Burlisch-Stoer integrator, to in- 
tegrate the regularized equations of motions. 

Computer codes for studying the dynamics of few body 
systems and star clusters or planetary systems are available. 
A very valuable r eview of th is kind of software industry is 
given in lAarsethl il999L 12003). But, interestingly, there is a 
lack of software for calculating the dynamical decay of sys- 
tems with a few < N < four dozen stars. Our endeavour 
is to fill this gap by a sophisticated software tool allowing 
us to efficiently study the decay of hierarchical and non- 
hierarchical configurations of some tens or hundreds stars 
down to the last remaining hard binaries or hierarchical 
higher-order multiple-stars, with the long-term-aim of em- 
bedding CATENA in a general-purpose A^-body code. 

An error analysis for the present application is provided 
in Seed 



The probability to draw 20 (P(fc < 20)), 10 (P(fc < 10)) or fewer 
stars heavier than 5 Mq from a Kroupa-IMF, the expectation 
value fi of the number of stars heavier than 5 Mq and the to- 
tal number of stars ntot (equivalent to the number of repeated 
experiments) are calculated for two different total cluster masses 
and two different physically possible upper stellar mass limits. 



ent ONC-masses and two different physically possible upper 
stellar mass limits is given in Tab.|!|] It is extremely unlikely 
that only ten stars have formed in the ONC if the IMF is 
universal. 

We note that the same argument can be applied to 
a more-evolved population: In an exploration of the full 
stellar populati o n of the Upper Scorpius OB association, 
IPreibisch et al.l i2002T) determined a total stellar mass of 
2060 Mq covering a volume of 35 pc in diameter. For the su- 
pernova progenitor they deduced a mass of m 40—60 Mq . An 
IMF steeper than 2.3 in the regime of massive stars would 
not have lead to the formation of such a massive star in the 
young star-cluster-stage of the Upper S corpius OB associa- 
tion 5 Myr ago for this mass of 2060 Mq iWeidner fc Kroupal 
2006). This further supports that the IMF may not be 
steeper than 2.3 for massive stars. IPreibisch et alJ ll2002l) 
listed 19 stars heavier than 5 Mq. This is approximately 
half of the expected number of formed stars more massive 
than 5 Mq and constitutes the same problem as for the ONC 
due to similar initial cluster masses. Therefore, it can be ar- 
gued that O and B stars may have been ejected from their 
star forming region very early after their formation. 

So two questions arise assuming the IMF is invariant: 
Why does the Trapezium still exist and where are the miss- 
ing OB-stars? 



4 INITIAL CONDITIONS 

To address the questions mentioned above we investigate 
three models which consist of the stars listed in Tab.0 

In the first model we study the stability of the actually 
observed Trapezium system consisting of O 1 A, O 1 B, O C 
and O 1 D precisely. In the second model, it is assumed that 
all currently observed OB stars in the ONC (Tab. [SJ were 
initially in a compact configuration as a core at the centre 
of the ONC, due either to mass segregation or ab-initio. In 
the third model we start with an OB core coming close to 
the expected number of 38. To find a suitable set of stars, 
all presently observed OB stars are used four times giving 
40 stars (4 times Tab. |SJ. 

The compact settings of OB-stars are motivated by the 
outcom e of the analytical investigation by Clar ke~fc Pringlel 
jl992j) that ma s sive stars form in compact groups. 
iBonnell fc Daviesl (|l998) concluded that the positions of 
massive stars in the Trapeziums cluster in Orion cannot be 
due to dynamical mass segregation, but must have formed 
in, or near, the centre of the cluster. 

For each of these three models 1000 configurations are 
created where the stars from Tab. are uniformly dis- 
tributed over a sphere with t he compact Trapezium radius of 
0.025 pc <HiUenbrandll997l) . The velocities are drawn from 
a Gaussian distribution with a velocity dispersion resulting 
from the virial theorem, 

G A/tot. ob 

a = R (6) 

(Tab. 01- After this the velocities are re-scaled slightly to 
ensure initial virialisation. 

This simple model does not include the rest of the ONC. 
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To estimate its effect on the core decay the ratio of the 
internal and external forces can be calculated. The OB-star 
core of radius r consists of n stars having the mean mass m. 
The gravitational force on one star is then 



Table 4. Force ratio for the 4-, 10- and 40-body model 



F n = G- 



(7) 



The force exerted by the rest of the ONC on one star in the 
core can be estimated by the Plummer force 



F p i = GmM p i(r 2 + b 2 )-ir 



(8) 



assuming the cluster can be represented reasonably well by 
Plummer model, which has been shown to be the case 



The resulting force ratio is 



$ = 



F P i 



nm 
M p i 



1 + 



(9) 



The mass M p \ is the cluster mass minus the mass of the 
OB-stars. The resulting force ratios can be seen in Tab. 0] 
The core dynamics is dominated by its self-gravitation. 

The escape velocities for the isolated core and the total 
Plummer sphere can also be compared. Both are obtained 
from the conservation-of-energy-theorem. The escape speed 
from the centre of the Plummer sphere, w e , P i, and the escape 
speed from the surface of an isolated OB-core are given by 



fe.pl — 



2Gm c 



«e,OB = 



2GmoB 



(10) 



respectively, where ro — 0.025 pc is the initial radius of the 
OB-core. For the 4- and 10-body model the escape speeds 
for the isolated model and the true embedded situation are 
comparable. In the 40-body model the escape speed from 
the core is dominated by the core itself. 

A second issue associated with the cluster shell of low- 
mass stars is two-body relaxation between an OB-star and 
the low mass stars of the cluster. Energy may be transfered 
from the OB-star core and ejected or evaporated OB-stars 
to the rest of the cluster. The relaxation time of the ONC 
is about 18 Myr llKroupall2005l) . The relaxation time for a 
heavy star is given by multiplying the relaxation time with 
the ratio of t he mass of th e most massive star and the mean 
stellar mass jSpitzerlll987Tl and describes the time-scale of a 
massive star to sink towards the cluster centre, 



irelax.OB 



moB 



(11) 



where the average mass m of a star is 0.35 Mq using 
a Kroupa-IMF. The resulting energy transfer time-scale 
ranges from 0.14 Myr (45.7 M s ) up to 1.26 Myr (5 M ), 
thus being shorter or comparable to the time spanned by 
the simulations and therewith probably an important issue 
in our context, given the age of the ONC ~ 1 Myr. In the 
case of no equipartition instability, energy transfer stops af- 
ter reaching energy equipartition, 

m < v 2 >= moB < Vq B > , (12) 

where < v 2 > (< Vq B >) is the mean square velocity of the 
mean-mass stars (OB-st ars, respectively). U sing a velocity 
dispersion of 2 km s" - 1 faillenbrandl 1 1997ft for the mean- 
mass stars, the relation above and the energy theorem it 
can be calculated that the velocity of a 5 Mq (45 Mq) 
is low enough such that the movement of the OB-stars is 



n 


Af OB /M 


m/M 




«e,OB 


Mcl/M© 


v e,pl 


4 


88.4 


22.1 


94.5 


5.6 


1721.4 


7.1 


10 


167.2 


16.7 


178.6 


7.7 


1800.2 


7.3 


10 


668.8 


16.7 


714.2 


15.4 


2301.8 


8.2 



n is the number of stars the model consists of, Mob (cf. Tab.0 is 
the total mass contained in the OB-stars, m is the mean mass of 
an OB-star, $ is the resulting force ratio using a Plummer mass 
of 1633 Mq . Given the observed core ra dius of the ONC of about 
0.19 pc iHill cnbrand & Hartmann 199^) the related Plummer pa- 
rameter of the ONC is about 0.3 pc. ti c ,OB is the escape speed 
in km/s from the surface of an OB-core with radius of 0.025 pc, 
i> e , p l is the escape speed in km/s from the centre of a Plummer 
sphere with mass M c \ = 1633 Mq + Mob- 



Table 5. Identity of the stars used in the three models (Tab. IT] 



Name 



Parenago SpT m/Mg 





1865 


09V 


18.9 


q 1 b 


1863 


B0V 


7.2 


e 1 c 


1891 


07V 


45.7 


e x D 


1889 


BOVp 


16.6 


o 2 a 


1993 


09V 


31.2 


o 2 b 


2031 


B1V 


12.0 


LP Ori 


1772 


B2V 


7.2 




1956 


B3 


6.1 


NU Ori 


2074 


B1V 


16.3 


HD37115 


2271 


B5V 


5.7 



Stellar data for all OB-st ars over 5 Mq gi 
Spectral type after Ivan Altena et alJ 



ivcn by Hillenbrand 
1988 



constrained to be within a radius of 0.026 pc (0.0084 pc). 
So the current observed OB core has an extension con- 
sistent with energy equipartition. Following iHeggie fc Hud 
( 2003) the heavy stars are so concentrated that the lighter 
stars have been expelled from the core and they no longer 
have a significant role. This is also suggested by the ob- 
served deficit of low-mass star s in the core of the ONC 
jHillenbrand fc Hartmannlll998lb 

We conclude that the effect of two-body relaxation be- 
tween low-mass stars and the OB-stars may be of minor im- 
portance and that these simulations suffice to demonstrate 
the time-scale problem of the ONC, and that the OB-star 
core-decay-model may explain the OB-star number problem 
of the ONC. While full-scale iV-body calculations capture 
the entire relevant physics, our approximations allow us to 
compute a very large number of renditions (here 5000 in to- 
tal) which is necessary given the low frequency of massive 
stars. Future A r -body calculations of individual set-ups will 
be used to check our results. 



5 FINDING TRAPEZIUM SYSTEMS 

We define a trapezium system to consist of a few stars having 
pairwise distances of the same order. Here the whole system 
is scanned to determine the maximum number of stars in 
a configuration in which the pairwise distances lie between 
two boundaries: When studying the stability of the ONC 
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a h 
b ^ 



Figure 1. Illustration of the trapezium-finding routine (see Sec, 151 
for explanation). 

Trapezium these boundaries are 0.01 and 0.05 pc. When 
studying the total OB-star distribution these boundaries are 
and 0.05 pc. 

This procedure is illustrated in Fig. Consider a con- 
figuration consisting of six bodies. A table containing the 
pairwise distances is created. All subsets of particles having 
a pairwise distance between two boundaries are determined. 
Of all these subsets the one having the most members is 
the extracted trapezium system. For the explicit example 
above, assume that a trapezium system of dimension a is 
searched. The pairwise distances may be allowed to devi- 
ate by 20 per cent from this dimension. Then the subsets 
found are: {2, 3}, {2, 4}, {2, 5}, {3, 4}, {3, 5}, {4, 5}, {2, 3, 4}, 
{2, 3, 5}, {2, 4, 5}, {3, 4, 5} , and {2, 3, 4, 5}. So the trapezium 
system consists of four bodies. 

If a trapezium system of the dimension b is searched 
all subsets extracted from the distance table are: {1,2,6}, 
{1, 3, 6}, {1, 4, 6}, and {1, 5, 6}. Four candidates for a trapez- 
ium system of dimension b are found. But they all have the 
same number of members. 

Based on this algorithm a set of bodies can contain 
no trapezium system of a certain dimension, or a trapezium 
system can have two or more members. But it is not possible 
to find a trapezium system consisting only of one body. 

If the number of stars within a certain sphere is of in- 
terest the lower boundary must only be set to zero. 

The pairwise distances of the Trapezium stars in the 
ONC Ol lie approximately between 0.02 and 0.05 pc, 
whereby all of the OB stars can be fo und within nearly 1 pc 
radius around 01 toillenbrandlll997l) . 



6 DECAY OF OB-STAR CORES 

All 3000 configurations are integrated over 2 Myr. 95% of 
all runs have a relative energy error lower than 10 -14 in the 
case of the 4-body model, 10 -12 for the 10-body model and 
10" 10 for the 40-body model. 

6.1 Four-body model 

In the top-diagram of Fig. the decay curve of a four- 
body trapezium consisting of the four stars of O 1 Ori is 
plotted. After 1 Myr only 2.5 stars on average are found 
with a pairwise distance less than 0.05 pc, and 1 star on 
average is found with a pairwise distance between 0.05 and 
0.01 pc, where in both cases the 3d curves lie slightly below 



the 2d-projection decay curve. This demonstrates the time- 
scale problem pointed out above for the observed Trapezium 
which is marked by the bold x. It can be argued that these 
are only mean values which are gained from a number dis- 
tribution, and that a compact trapezium as is observed can 
survive for 1 Myr with some probability. That this is not 
the case can be seen in Fig. [3] where the distribution of the 
member numbers is plotted in a histogram. In 62.2 % of all 
runs no stellar configuration is found with a pairwise dis- 
tance between 0.01 and 0.05 pc, while in only 0.4 % of all 
runs a trapezium of the observed size is found after 1 Myr. 

On a first view this stands in contradiction to 
lAllen fc Povedal lll974f) . which is the only known work about 
stability of trapezia systems. Their result was that 63 % 
of all trapezia remain as a trapezia after 30 crossing times 
(~ 1 Myr). The reason is that they used a slightly different 
definition of a trapezium system: "let a multiple star system 
(of 3 or more stars). . . if three or more such distances are of 
the same order of magnitude, then the multiple system is of 
trapezium type.. . .Two distances are of the same order of 
magnitude, in this context, if their ratio is greater than 1/3 
but less than 3". A four-body system satisfying our defini- 
tion of a trapezium system, i. e. all pairwise distances are 
of the same order of magnitude, can evolve into a system of 
two singe stars and one close binary. Then three distances 
are of the same order of magnitude, i. e. the distance be- 
tween the two single stars, and the distances between each 
single star and one component of the binary. Their definition 
will detect a trapezium system. Our algorithm also detects 
a trapezium system but consisting only of three stars and 
having a different size than searched for. So the criterion 
by Allen and Poveda does not take the number of members 
of the trapezium system into account nor the size of the 
trapezium. For a quantification of the stability of a four- 
body system, the crucial point is the number of stars the 
trapezium consists of. The Allen-and-Poveda trapezia con- 
sist initially of six stars, and only one of all 30 configurations 
(3.3 %) retains the initial size after 1 Myr, and consists fi- 
nally only of four stars. Indeed in all their 30 runs binaries 
form consisting preferentially of the two most massive stars, 
therewith being quite consistent with our result. 

6.2 10-body model 

As in the four-body model the mean number of stars within 
0.05 and 1 pc is plotted in Fig. [5] The decay of an initial 
ten-body OB star core can neither reproduce a four-body 
trapezium with a diameter of 0.05 pc nor the entire present- 
day ONC OB population (Fig.|3J. 

6.3 40-body model 

If it is assumed that the ONC had an initial OB-star con- 
tent of nearly forty as expected from the canonical IMF 
combined with the estimated stellar mass of the ONC of 
about 2200 Mq then the remaining number of OB stars af- 
ter 1 Myr within a sphere of 1 pc radius comes close to the 
observed value of ten (Fig. |5J. The probability to observe a 
trapezium at an age of 0.5 Myr is 13.8 % but only 0.7 % at 
an age of 1 Myr (Fig.[3J. But counted together with the sys- 
tems containing more than four stars the probability to find 
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a compact trapezium increases to 0.7 + 2.9 = 3.6%. Note 
that the estimated ONC mass is comprised of the observed 
mass (1800 Mq) plus the estimated mass in missing 30 OB 
stars (400 Mq). The true initial ON C mass may have been 
about 4000 M llKroupa et al.ll200ll) . 

In Fig. 2]the spatially cumulative and velocity distribu- 
tions of the O and B stars after 1 and 2 Myr are plotted. 
After 1 Myr more than 75 % (30 of 40) of all stars have 
larger distances to their common centre of mass than 2 pc, 
and after 2 Myr 75 % of all stars are more than 4 pc away 
from their centre of mass. Therefore only ten of forty stars 
heavier than 5 Mq remain at the cluster centre as observed 
in the ONC. It can be seen that more O stars than B stars 
are at very large distances as well as more O stars than B 
stars have very high velocities. This comes about because 
initially B stars tend to evaporate by energy redistribution 
rather than being ejected by close encounters. Because O 
stars are heavier they form tighter configurations than B 
stars and they are then involved in close binary interactions 
le ading to high ejec t ion ve locities. This confirms the result 
of IClarke fc Pringlj lll992l) qualitatively. 

Given the spatial distribution of OB-stars after an evo- 
lution of 1 and 2 Myr, an extrapolation to an evolution age 
of 5 Myr predicts that nearly 50 % of all OB-stars cover a 
volume of 25 pc in diameter which comes close to the ob- 
served properties of the Upper Scorpius OB association as 
pointed out in Sec. 12.21 



7 ERROR ANALYSIS 

An important question concerning iV-body simulations is 
whether we can trust them. This question arises from the 
fact that the basic process for the decay of iV-body systems 
are close and highly eccentric encounters of stars with a 
consequent redistribution of energy. These kind of encoun- 
ters are the most important sources for orbit-errors. Due 
to the exponential instability the nu merical and true solu- 
tion deviate increas ingly with time. IIGoodman et al.lll993t 
iHeggie fc Hutll2003l) 

As there exists no general analytic solution for systems 
with more than 2 bodies the integrals of motion have to 
be checked for conservation to control the integration. The 
most sensitive quantity is the total energy. Smaller step sizes 
reduce the amount of the error but prolong the duration of 
the integration, while step-sizes that are too small lead to 
the accumulation of roundoff errors. So a balance between 
effic iency and accuracy need s to be found. 

iDeionghe fc Hutl ll 198d) determined the accuracy of the 
integration of 14 small-iV configurations using the time- 
reversal-test co mpared with en ergy conservation, leading to 
the conclusion <IAarsetbl l2003) that using the energy error 
only is questionable for establishing exact integrations. 

When using the statistical approach, high accuracy is 
not essential to obtain meaningful resu l ts, provided the s am- 
ple is sufficiently large jAarsethll2o"o3T) . IValtonenl dl974h de- 
termined that the distributions of eccentricity, terminal es- 
cape velocity and life-time in 200 3D experiments did not 
show any clear accuracy dependence for a relative energy 
error range from 5 • 10~ 4 to 3 ■ 10~ 2 . 

Therefore, despite the fact that the orbits are com- 
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Figure 2. Decay curves for the N =4, 10,40 models. Plotted is the 
mean number of stars for configurations with maximum member 
numbers having a certain pairwise distance as a function of time. 
The bold x marks the position of the present ONC-TS with an as- 
sumed age of 1 Myr, and the bold o marks the observed OB-stars 
in the ONC. top: Decay of the 4-body model where the pairwise 
distance lies (from bottom to top) between 0.01 and 0.05 pc (3d), 
between 0.01 and 0.05 pc (2d), below 0.05 pc (3d) and below 
0.05 pc (2d), middle: Decay of the 10-body model where the pair- 
wise distance lies (from bottom to top) below 0.05 pc (3d), below 
0.05 pc (2d), below 1 pc (3d) and below 1 pc (2d), bottom: Decay 
of the forty-body model where the pairwise distance lies (from 
bottom to top) below 0.05 pc (3d), below 0.05 pc (2d), below 
1 pc (3d) and below 1 pc (2d). 
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Figure 3. Distribution of the remaining number of stars in 
trapezium configurations having a pairwise stellar distance be- 
tween 0.01 and 0.05 pc in 2-d projection after 0.5 and 1 Myr. top: 
4-body model, middle: 10-body model, bottom: 40-body model. 



pletely wrong after many crossing times, the statistical out- 
come from many equivalent A^-body experiments is reliable. 

To test what value of energy error is acceptable we run 
all 1000 four-body configurations three times with different 
step-size parameters. The resulting mean energy errors are 
5.65 ■ 10~ 12 , 2.73 • 10~ 6 and 6.19 ■ 10~ 2 with increasing step- 
size parameter (Fig. |SJ . 

To compare the statistical error with the numerical 
error we interpret this analysis as a series of Bernoulli- 
experiments. One experiment is the determination of a 
trapezium consisting of A/stars stars after a certain time T. 
The outcome can be yes with probability p and no with the 
probability 1 — p. This experiment is repeated n = 1000 
times. Therefore the probability to get fc times the event yes 
is given by 



w - G) 



k /i \n — k 

P (1-p) 



(13) 



The event probability is approximated by the mean proba- 
bility 



p = p = (fci + k' 2 + fc 3 /3)/n = k/n, 



(14) 



where fci is the number of events in each of the three ex- 
periment series. The variance of the Bernoulli-distribution 
is given by 




0.001 



30 ,, ,40 
v/ km/s 



Figure 4. top: Cumulative spatial distribution of O and B stars 
for the 40-body model after 1 and 2 Myr measured relative to 
the total centre of mass, bottom: logarithmic Distribution of the 
velocities of O and B-stars for the 40-body model after 1 Myr. 
Double stars are considered using their centre of mass velocity. 



Afc = np - 



(15) 



The corresponding probability p, number of experiments n, 
mean number of events fc, variance y/ Afc 2 and one sigma 
errors A% for the histograms in Fig. are given in Tab. [(J 
As an example we consider trapezia consisting of 2 bodies 
after 1 Myr. The mean number of runs having a trapezium 
of 2 bodies after 1 Myr is 802 out of 1000 (80.2 %). The one 
sigma variance is 2.8 %. So all three experiments (81.9 %, 
79.9 % and 78.9 %) lie within 1 sigma around the mean value 
(80.2 % ± 2.8 % = 77.4 % - 80.2 %). So if a numerical error 
arises from the different choice of the step-size parameter it 
is not larger than the statistical error. We conclude that we 
can trust these A^-body simulations. 
From the virial theorem 



E = V/2, 



(16) 



where E is the total energy and V is the potential energy, 
the relative energy error is 



AE/E ss AV/V, 
and with 

M 2 
R ' 



V 



-G 



(17) 
(18) 
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step size: CTRL = 10 in CATENA 
mean energy error: 5.65 ■ 10 -12 
after 1 Myr: 
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step size: CTRL = 500 in CATENA 
mean energy error: 6.19 • 10~ 2 
after 1 Myr: 
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Figure 5. Error analysis using the example of the four-body 
decay determining the maximum number of members in configu- 
rations with a pairwise distance below 0.05 pc in 3d after 1 Myr 
(left column) and 2 Myr (right column) for three different step 
size parameters resulting in mean energy errors of 5.65 ■ 10~ 12 
(top), 2.73 ■ 10 -6 (middle) and 6.19 ■ 10 -2 (bottom diagram). 



Table 6. One sigma errors for the four-body decay (Fig. |SJ in- 
terpreted as Bernoulli-experiment. 
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Resulting variances for the error histograms in Fig.[^]if the search 
for a trapezium consisting of Altars is interpreted as a Bernoulli- 
experiment. 



AE/E w AR/R, 



(19) 



follows. This means that uncertainty in the pairwise dis- 
tances is of the order of the mean energy error. As the en- 
ergy errors are smaller than 6 • 10 -2 the distance errors have 
only a very slight effect on the number statistics. 



8 CONCLUSIONS 

We have pointed out that in the case of the Orion Nebula 
cluster two main problems exist: Its short decay-time im- 
plying the question as to why it exists, and the significant 
number of missing OB stars implying either that they have 
been lost if the IMF is invariant, or that the IMF had a 



highly non-standard Q3 > 2.7 which is unlikely because no 
other stellar population in a cluster with such an IMF is 
known to exist and the observed most-massive star in the 
ONC is significantly larger than that expected for this non- 
standard IMF. It is extremely unlikely (3 out of 1000 cases) 
that a compact Trapezium system consisting of four stars 
can survive for more than 1 Myr. The assumption that the 
initial number of OB stars was about 40 increases the prob- 
ability to observe a Trapezium system after 1 Myr (36 out 
of 1000 cases). We infer that the ONC Trapezium system 
could be an OB-star core in its final stage of decay. 

This scenario is supported by the fact that the spa- 
tial distribution of these forty OB-stars obtained from the 
numerical simulations after 5 Myr comes close to the ob- 
served spatial distribution of the OB-stars in the Upper 
Scorpius OB association which is assumed, using the total 
mass, to have had the same young star-cluster progenitor as 
the ONC. 

Given its total mass the ONC indeed ought to have four 
times as many OB stars than are observed, namely 40. This 
suggests that about 30 OB stars may have been expelled 
from the ONC if the IMF is invariant. Starting with an 
initial number of OB stars of about 40, stars are ejected 
due to three-body encounters, so that this model matches 
the observed number of OB stars in the ONC and the ONC 
Trapezium. 

As it has been shown that the probability that only 10 
OB-stars have formed in the ONC is very small, the missing 
stars must be somewhere. After 1 (2) Myr 90 (70) % of all 
OB-stars should be found within a 10 pc radius around the 
ONC centre. It is straight-forward to search for these missing 
OB-stars in OB catalogues. If they cannot be found then the 
IMF may be steeper at the high mass end. Or OB stars form 
with an initially high binary fraction so that ejections occur 
faster and the resulting velocities are higher, placing them 
even further away from the ONC centre after 1 (2) Myr. 

Due to the absence of primordial binaries ours is a 
conservative result because binaries enhance the ejection 
rates. Therefore the influence of primordial binaries must 
be investigated in further experiments, which will also 
need to consider gas expulsion from the e mbedded cluster 
jKroupa et alj|200ll: IVine fc Bonnelill2003l) . IVine fc BorinelJ 
(2003) investigated the evolution of cores of young star clus- 
ters and their massive stars but used a smoothed potential 
to avoid the difficulties coming up with the occurrence of 
close encounters. But these close encounters are the energy 
sources for massive star ejections and the reasons why young 
stars can be found far away from the cluster centre within 
a short period of time. 

The influence of the cluster potential and especially of 
two-body relaxation with low-mass stars in the cluster shell 
may also be investigated in further simulations. Of what 
kind this influence is is still somewhat unclear. Stronger con- 
straining forces by the cluster potential and energy-loss by 
two-body relaxation may return some of the OB-stars evap- 
orated with low velocities from the core. This may stabilize 
the core on the one hand, but may also lead to a faster de- 
cay by shifting the velocity spectrum of the ejected stars 
to higher velocities because the stellar density at the centre 
would be higher. This would increase the probability of close 
encounters and high velocity ejections. If relaxation indeed 
stabilises the core then more OB-stars are expected to re- 



An abnormal massive-star- ONC-IMF? 9 



main in the ONC worsening the OB-stars discrepancy. Ours 
is therefore a conservative result. 

As a final note, Kroupa et al. (2001) have presented 
star-cluster-formation calculations that reproduce the ONC 
at an age of 1 Myr and the Pleiades at an age of 100 Myr. 
These models, however, are about 1.8-times as massive as 
the ONC mass used here (2200 M©) implying that if the IMF 
was canonical then the ONC may have had 1.8 x 40 = 72 
stars more massive than 5 M© . This would pose an increased 
challenge, because as is evident from Fig. 01 40 OB stars al- 
ready lead to an acceptable match with the data, so 72 would 
increase the probability of finding a trapezium-configuration 
at an age of 1 Myr, but would lead to too many OB stars 
within the cluster (as can be deduced from the lower-panel 
in Fig.HJ. 

Clearly, such models with a high initial multiplicity frac- 
tion need to be constructed for further studies of the intri- 
cate interrelation of the IMF with stellar dynamics in young 
clusters. 
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APPENDIX A: A PRACTICAL NUMERICAL 
FORMULATION OF THE IMF 

Al The general IMF 

When doing research using the IMF a multi-power-law, 



,mi ow < m < ma 
,uih < m < mo 
( — ) ,m < m < mi 

m H / V mo / \ mi / ' — — 



CM 



rap 



(A2) 



with exponents for its canonical form, 

ao = +0.30 , 0.01 < m/Af© < 0.08, 

ai = +1.30 , 0.08 < m/M e < 0.50, 

a 2 = +2.30 , 0.50 < m/M e < 1.00, 

a 3 = +2.35 , 1.00 < m/M© < +oo, 

jKroupa et alJ Il993l: iReid et alJ l2002t iKroupal l200ll: 
IWeidner fc Kroupal I2OO4F requires many if-statements in 
code implementations. Note however that the canonical IMF 
parametrisation constitutes a two-part power-law IMF in 
the stellar regime: «i = 1.3 for m < 0.5 M© and 02,3 = 2.3 
for m > 0.5 M©. Here we present a general formulation 
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for an IMF that avoids these difficulties with IF-statements 
and allows any number of mass segments and types of in- 
terpolating functions. Complicated code constructions con- 
taining if-statements are replaced by two easy loops. The 
described method is realised by some very handy functions 
implemented in a shared C-library, and is available on re- 
quest or at the AlfA homepage 1 . The formulation described 
below can be applied to any arbitrary distribution functions 
for any purpose. 

Due to historical reasons multi-power-law IMFs start 
indexing intervals and slopes at zero instead of one. For sim- 
plicity we here index n intervals from 1 up to n. 

Consider an arbitrary IMF with n intervals fixed by the 
mass array [mo, . . . , m n ] and the array of functions ft,. . . ,/ n . 
So on the i-th interval [mi_i, mi] the IMF is described by the 
function fi. For the case of a multi-power law the functions 
may be chosen to be 



Mm) 



(A3) 



To make this power-law description correspond to the multi- 
power-law given above the IMF-slope indices above need to 
be shifted by one. 

The following general IMF-description does not require 
power- laws for the functions fi, but also any kind of function 
is allowed. This incl udes log-normal IMFs jMiller fe Scald 
ll979l:IChabrierll2003ft . 

With the two 0-mappings 



©[](*) 
and 

6] [(a:) 



1 x > 
x < 



1 x > 
x < 



the function 

T[i] (m) = ©[ ] (m — mi_i)0[ j (mi — m) 



(A4) 



(A5) 



(A6) 



can be defined. It is unity on the interval [mi_i,mi] and 
zero otherwise. The complete IMF can be conveniently for- 
mulated by 



£(m) = k Yl A(m - mj) ^ T [n (m) ¥, /i(m) 



(A7) 



where is a normalisation constant and the array 
• . ,$n) is to ensure continuity at the interval bound- 
aries. They are defined recursively by 



*i = 1 , = 



/i-i(mj-i) 
/;(m;-i) 



(A8) 



For a given mass m the rpi makes all summands zero ex- 
cept the one in which m lies. Only on the inner interval- 
boundaries do both adjoined intervals give the same contri- 
bution to the total value. The product over 



A(x) = 



0.5 x = 
1 x^O 



(A9) 



halves the value due to this double counting at the interval- 
boundaries. In the case of n equals one (one single power 

1 http://www.astro.uni-bonn.de 



law), the empty product has, by convention, the value of 
unity. 

An arbitrary integral over the IMF is evaluated by 

6 fb i>a 

f (m) dm = £(m) dm- f (m) dm , (A10) 

J mg J mo 

where the primitive of the IMF is given by 

£(m) dm = k S Oj [(a — mi) / fi( m ) dm 

i=l ^mi_i 

+ fc^r [i] (a)* i / /i(m)dm. (All) 

i=l Jmi_i 

The expressions for the mass content, i.e. m £(m), and 
its primitive are easily obtained by multiplying the above 
expressions in the integrals by m. 



A2 The individual cluster IMF 

A2.1 Normalising the IMF 

The IMF denotes the number of stars per mass inter- 
val. Therefore the normalisation depends on the clus- 
ter mass. Here we follow the normalisation strategy by 
IWeidner fc Kroupal (120041) . This method requires two fur- 
ther masses. m max * is the maximum physically possible stel- 
lar mass and m max is the expected maximum stellar mass in 
a given cluster of mass M c \- With £ = k £fc(m) two equations 
defining m max * and m ma x result: 



M r , = k 



1 = k 



m £fc(m) dm, 



£k(m) dm. 



(A12) 



(A13) 



To solve these two equations for k and m max they can be 
divided by each other leading to an expression for the cluster 
mass as a function of m max 

M c i= / m^(m)dm/ / £ fc (m) dm. (A14) 

J mo J m max 

As a function of m max it is continuous and strictly mono- 
tonically increasing and its image is 72.>o- Therefore the 
existence of a solution m max is concluded by a connectiv- 
ity argument and the uniqueness follows from the strict 
monotony. This solution can be gained using any equation- 
solving method. 

Above m max no stars can be found and the IMF in a 
star cluster can be expressed by 

Cci(ra) = 8[ ](m max - m) {(m) . (A15) 



A2.2 Dicing stars - the mass-generating function 

When doing research on the IMF using Monte Carlo simu- 
lations or in setting-up star clusters for A^-body simulations 
a finite set of random masses distributed according to the 
IMF have to be diced. A random number X is drawn from a 
constant distribution and then transformed into a mass m. 
The mass segments transformed into the X-space are fixed 
by the array Ao,. . • ,A n defined by 
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Xi = / fd(m) dm. (A16) 

If P(X) denotes a constant distribution between and A n , 
both functions are related by 

/ fci(m') dm' = / P(X') AX' = X (A17) 

J niQ J 

for a uniform distribution -P(A"). The solution of this equa- 
tion for m is given by 

n-1 

j=i 

where F\ is the primitive of f\, F^ 1 is the primitive's inverse 
mapping and \F[ are mappings which are unity between Ai-i 
and Ai and zero otherwise. 



APPENDIX B: FINDING THE NUMBER OF 
EXPECTED OB-STARS IN THE ONC 

The observed total stellar mass of the ONC may be less 
than the initial one if OB-stars have been ejected. The total 
stellar mass M c \ is related to the IMF by 

/* m max 

M c i = Af< 5 + / m £(m) dm, (Bl) 

where M<5 is the observed mass in stars less massive than 
5 M Q . For a total mass for the ONC of 1800 - 3300 M the 
expected maximum mass m max lies in the range 50 - 63 M Q 
JWeidner fc Kroupall2006l) . If 5 M < m max , which is the 
case, the IMF can be normalised directly by 

M< 5 = k / m £fc(m) dm, (B2) 



where mo = 0.01 Mq is the opacity-limited minimum frag- 
mentation mass. The maximum stellar mass is then deter- 
mined in the second step solving 



1 = / £(m) dm, (B3) 

J m max 

whi ch means there exists one most massive star in the clus- 
ter JWeidner fc Kroupall2004 . The expected number of OB 
stars is then given by 

/>m max 

A^ob = / eel dm. (B4) 



